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Foreword 

This report has been prepared to expedite early dissemination of the information 
generated under the contract. The data and conclusions must be considered preliminary and 
subject to change as further progress is made on this program. This is a progress report covering 
the work done during the fourth 12 months of the contract; it is not a final report. The NASA 
Program Manager is Dr. C.C. Chamis. 
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1.0 Introduction 


This technical program is the work of the Life Analysis and Methods Technology Section 
of GE Aircraft Engines in response to NASA RFP 3-537260, “Coupled 
Structural/Thermal/Electromagnetic (CSTEM) Analysis/Tailoring of Graded Composite 
Structures.” The overall objective of this program is to develop and verify analysis and tailoring 
capability for graded composite engine structures taking into account the coupling constraints 
imposed by mechanical, thermal, acoustic, and electromagnetic loadings. 

The first problem that was attacked is the development of finite elements capable of 
accurately simulating the structural/thermal/electromagnetic response of graded composite 
engine structures. Because of the wide diversity of engine structures and the magnitudes of 
the imposed loadings, the analysis of these is very difficult and demanding when they are 
composed of isotropic, homogeneous materials. The added complexity of directional 
properties which can vary significantly through the thickness of the structures will challenge 
the state of the art in finite element analysis. We are applying AE’s 25 years of experience in 
developing and using structural analysis codes and the exceptional expertise of our University 
consultants toward the successful conclusion of this problem. To assist in this, we drew heavily 
on previously funded NASA programs. 

We built on NASA programs NAS3-23698, 3D Inelastic Analysis Methods for Hot Section 
Components, and NAS3-23687, Component Specific Modeling, in the development of the 
plate and shell elements. In addition to these two programs, we drew on NAS3-22767, Engine 
Structures Modeling Software System (ESMOSS), and NAS3-23272, Burner Liner 
Thermal/Structural Load Modeling, in Task III when we generated a total CSTEM Analysis 
System around these finite elements. This guarantees that we are using the latest computer 
software technology and produced an economical, flexible, easy to use system. 

In our development of a CSTEM tailoring system, we built on NASA Program 
NAS4-22525, Structural Tailoring of Engine Blades (STAEBL) and AE program. Automatic 
Improvement of Design (AID), in addition to the program system philosophy of ESMOSS. 
Because of the large number of significant parameters and design constraints, this tailoring 
system will be invaluable in promoting the use of graded composite structures. 

All during this program, we availed ourselves of the experience and advice of our Low 
Observables Technology Group. This will be particularly true in the Task V proof-of-concept. 
Their input will be used to assure the relevance of the total program. 

Figure 1 shows our program and major contributions in flowchart form. This gives a visual 
presentation to the synergism that will exist between this program and other activities. 

Figure 2 depicts an integrated analysis of composite structures currently under 
development in the composite users’ community. The severe limitations of such a system are 
not highlighted because three major steps in the process are not shown. Figure 3 adds these 
steps. The analysis system really begins with a definition of geometry. A user then defines a 
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Figure 1. Program Flowchart. 
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Figure 2. Composite Analysis System. 
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Figure 3. Total Composite Analysis System. 


finite element model simulating this geometry and the anticipated loading. The process then 
moves to defined Step 3. One cycle through the process ends with the prediction of individual 
ply average stresses and strains. Now comes a significant productivity drain, namely, manual 
intervention to evaluate these stresses and strains against strength and durability limits. Based 
on this, the user must decide to (1) change the finite element model, (2) change the composite 
laminate, (3) both of the above, or (4) stop here. 

Obviously, there is a considerable cost savings to be obtained by selecting Number 4. The 
CSTEM system will obviate the reasons for selecting Number 4. This system, shown in Figure 
4, begins with the definition of geometry, as before, but then proceeds to a definition of master 
regions which contain all of the necessary information about geometry, loading, and material 
properties. Step 3 is a constitutive model which develops the necessary structural, thermal, and 
electromagnetic properties based on a micromechanics approach. Furthermore, this 
constitutive model contains the logic to generate the global finite element model based on the 
variation of the properties, as depicted in Figure 5. Using a nonlinear incremental technique, 
those global models are solved for their structural, thermal, and electromagnetic response. 
Based on this response the global characteristics are evaluated, with convergence criteria and 
decisions made on remodeling. Once the global characteristics meet the accuracy 
requirements, the local characteristics are interrogated and decisions made on remodeling 
because of strength, durability, or hereditary effects. Once this cycle has been stabilized, 
optimization is performed based on design constraint. 

1 .1 Executive Summary 

“CSTEM” is the acronym for the computer program being developed under the NASA 
contract, “Coupled Structural/Thermal/Electromagnetic Analysis/Tailoring of Graded 
composite Structures.” The technical objectives for this program are to p'Muce radar signal 
transparent structures having high structural performance and low cost. The multidisciplines 
involved are all highly nonlinear. They include anisotropic, large deformation structural 
analysis, anisotropic thermal analysis, anisotropic electromagnetic analysis, acoustics, and 
coupled discipline tailoring. The CSTEM system is a computerized multidiscipline simulation 
specialized to the design problems of radar absorbing structures. The enabling technical 
capabilities are implemented in a special 3D finite element formulated to simultaneously tailor 
the geometrical, material, loading and environment complexities of radar transparent 
structures for cost effective optimum performance. 

In each enabling technical discipline a decoupled stand-alone 3D finite element code has 
been developed. An executive program with controlling iterative solution techniques performs 
the nonlinear coupling among the participating technical disciplines. A geometry and finite 
element model generator specialized for graded composites has been developed as an intimate 
part of this analysis system. 

The structural analyzer is built around the 8-, 16-, and 20-noded isoparametric finite 
elements with emphasis on the 20 noded. A graded composite constitutive model has been 
developed for these elements which uses the NASA Lewis program ICAN as a subroutine to 
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Integrated analysis of composite structures 
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Figure 5. Constitutive Model • Structural Model Interaction. 



perform composite micromechanics and supply, to the analyzer, the requisite composite 
properties. The composite stiffness gradient controls the finite element definition of a 
structure with two major parameters to vary the number of elements through the thickness 
and the number of numerical quadrature points within an element. A unique set of local 
stiffness characteristics is developed for each numerical integration point. Integration of these 
local characteristics over the volume of the element provides total element simulation of 
composite structures including such effects as twist-bend coupling. 

The structural analyzer also performs large deformation analysis using a unique 
incremental updated Lagrangian approach with iterative refinement. Testing of this capability 
against classical large deformation problems has shown it to be both more accurate and more 
economical than available alternatives. Connected with this technical capability is a deformed 
position eigen-analyses capability. All or selected portions of the nonlinear stiffness terms can 
be incorporated into these eigenanalyses. This capability has been checked out against 
available test data and other computer codes. 

These capabilities have been combined with an optimizer to perform a totally automated 
integrated analysis of graded composite structures. A final task in the program will be a design 
demonstration of the tailoring capabilities of the CSTEM system. 

In order to reach the thermal analysis goals set for CSTEM, the same 8-, 16-, and 20-noded 
isoparametric finite elements are utilized. Four heat transfer solution options are available: 
linear steady state, nonlinear steady state, linear transient, and nonlineaar transient. To 
overcome the previous economic penalty associated with finite element vis-a’-vis finite 
difference heat transfer, a unique solution technique is employed. This is a Newton-Raphson 
iterative technique with right hand side pseudo-fluxes. The code will perform the heat transfer 
analysis of a thermally anisotropic material considering conduction, convection, and radiation. 
Table 1 lists the parameters involved. 

Routines for the calculation of absorption of electromagnetic waves have been written and 
checked out. The first method developed was based on a data bank of absorptivity values for 
given material types specified at discrete values of temperatures, frequency, and polarization 
angle. The information for a specific material or materials is read from the data bank file into 
arrays which are used by the program. The absorptivity is linearly interpolated from these 
discrete values to the local values of temperature, frequency and polarization angle. 

Calculations are made for one given frequency, orientation, and path of an 
electromagnetic wave at a time. Multiple frequencies and/or different paths travelled through 
the structure are handled by separate calculations for the different parameters. 

A wave coordinate system is associated with the electromagnetic wave. This wave 
coordinate system is defined such that the direction of propagation is along the positive y axis 
and polarization is measured from the positive x axis. 
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Table 1. Thermal Analyzer. 
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The element face upon which the wave impinges is input and the angle of incidence of the 
wave is calculated from the dot product of the inward normal of the layer subsurfaces and the 
positive y axis of the wave coordinate system. The angle of incidence must be between 0 and 
90 degrees in absolute value. The temperature is calculated using the element shape functions. 

The polarization angle is calculated using the dot product between the projection of the 
wave polarization on the layer midsurface and the material principle direction. The 
polarization angle is calculated as being between 0 and 90 degrees. 

As a second alternative, a method was developed that calculates the electromagnetic 
response and determines the reflected amounts of electromagnetic energy as well as the 
attenuation of the transmitted amount. This technique is based on Maxwell s equations, Snell s 
law, and the Fresnel formula. 

A third alternative has been provided by installing the computer program WAVES as a 
subroutine in CSTEM. This program calculates the reflection and transmission of 
electromagnetic waves given a stackup sequence of materials and their electromagnetic 
properties. 

The capability to determine acoustic characteristics due to structural vibration has been 
incorporated into CSTEM. The approach used determines the radiation efficiencies of a 
structure for each vibration mode as a function of frequency. An eigenanalysis produces the 
fundamental modes and mode shapes. Once the radiation efficiencies for each mode are 
calculated, the total sound power is obtained by a modal summation of the contribution from 
each mode. 

CSTEM tailoring capability has been built on the STAEBL (Structural Tailoring of Engine 
Blades) computer program obtained from NASA Lewis. This program consists of two major 
modules; CONMIN, which performs optimization and ANALIZ which supplies the parameter 
to be optimized. The CONMIN module was abstracted from STAEBL and coupled with the 
CSTEM structural, thermal, electromagnetic and acoustic application modules. This was 
successfully debugged, verified, and validated. 

Further improvements have been made to the CSTEM micromechanics routines and a 
unique modeling and solution technique for ply drop-off has been developed. 
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2.0 Technical Progress 


2.1 Demonstration of Coupled Solution Capability 

A simulated composite duct was generated to use in testing the coupled 
structural/thermal/electromagnetic solution capabilities of CSTEM. The CSTEM model 
generator was used to automatically develop the finite element model shown in Figure 6 . 
Thirty-two 20-noded graded composite finite elements were generated with four composite 
layers ( + 60°, -60° + 60°, -60°) in each element. The thermal problem consisted of convection 
at the inner surface from a temperature varying fluid and the outer wall being held at a constant 
temperature. Figure 7 plots the time variation of the convection coefficient and Figure 8 that 
of the internal environment temperature. The necessary thermal properties of the composites 
came from the internal ICAN data bank. The structural loading, in addition to temperature, 
consisted of internal pressure and an imposed duct end displacement. Figure 9 is a time history 
of the end displacement and Figure 10 is a time history of the internal pressure. For the 
electromagnetic absorption problem, simulated temperature, frequency, and polarization 
dependent absorption tables were supplied to the data bank (Figures 11 and 12). Figure 13 is 
a plot of the original and the loaded, deformed geometry. Figure 14 is a plot of the variation 
of the temperature through the duct thickness. Figure 15 plots the effective stress for a line of 
elements covering the distance from duct tip to duct root and showing the composite layer 
effect. Figures 16 and 17 show the variation around the circumference of the composite layer 
effective stresses and strain for the tip elements. Figures 18 and 19 show tlu same information 
for the root elements. The effect of the ± 60° layers is evident. This run required approximately 
89 CPU seconds on the CRAY-XMP. 

2.2 Alternate CSTEM Electromagnetic Technology 

A different method for calculation of electromagnetic response has been included in 
CSTEM that determines reflected amounts of electromagnetic energy as well as calculating 
the attenuation of the transmitted amount. Figure 20 shows an electromagnetic wave 
propagating in the + Z direction of the coordinate system attached to the wave. The electric 
field vector points in the + X direction and the magnetic field vector points in the + Y direction 
for a polarization of zero degrees. 

The wave equations which apply to these vectors can be obtained using Maxwell’s 
equations, and can be written as: 

v 2 H = 7 2 H and V 2 E = 7 2 E 


where the propagation constant, 7 , can be written as: 
y = a + jp 
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Figures. Composite Duct Model. 
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Figure 8. Environment Temperature. 
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Figure 9. End Displacement. 



Tim*, Mcond 


Figure 10. Internal Pressure. 
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Figure 12. Absorption at 600°F. 
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Figure 13. Orlginaland Loaded, Deformed Geometry. 
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Figure 15. Effective Stress Versus Distance from Duct Root. 
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Figure 16. Effective Stress at Tip Element Centroids. 



Figure 17. Effective Strain at Tip Element Centroids. 
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Figure 19. Effective Strain at Root Element Centroids. 
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Figure 20. Phase Relationship of Electromagnetic Wave. 


20 


Here, a is the attenuation factor and is found as: 


and p is the phase shift constant, written as: 

•’-“Vi (iMif.i) 

In the above equations, o> is the frequency of the electromagnetic wave in radians/second, 
is the permeability of the material in henries/inch, e is the permittivity of the material in 
farads/inch, and a is the conductivity of the material in siemens/inch. 

The index of refraction of a material is the ratio of the wavelength or phase velocity in free 
space to that in the material. Since the wavelength can be written in terms of the phase shift 
constant as: 



the index of refraction can be written in terms of the phase shift constant as: 



where p is the phase shift factor for the dielectric material and p 0 is the value of the phase 
shift constant for free space found as 

P = 6Vm £ 

o o o 

The calculation of the reflected and transmitted electromagnetic energies is done using 
Snell’s law of refraction and the Fresnel formulas. Snell’s law can be stated as: 


nisine = n 2 sin<|> 

where m and ni are the index of refraction of material 1 and 2, e is the angle of incidence of 
an electromagnetic wave propagating from material 1 to material 2, and $ is the angle of 
refraction of the transmitted portion of the electromagnetic wave. Figure 21 shows the 
geometry associated with an impinging electromagnetic wave. 

The Fresnel formulas assume an impinging wave with components of the electric vector 
that are in and out of the plane of incidence defined in Figure 21. This impinging wave is written 
as: 
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Figure 21 . Geometry of Snell's Law. 
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Ai = A p i (cos e I - sin e K) + Asi J 

where Api is the in plane amplitude and A s i is the out of plane amplitude, and e is the angle 
of incidence on the interface between materials 1 and 2 . The reflected portion of this wave can 
then be written: 

A2 = A p 2 (cos 0 1 + sin e K) + As 2 J 
where 

A - -A tan (8 - <t>) 

p2 pi tan (0 + <)>) 


_ sin (6 - d>) 

s2 si sin (0 + <J>) 


The refracted wave can be similarly written as: 
A3 = A p 3 (cos 4> I - sin 4> K) + A $3 J 
where 

. _ . 2 sin d> cos 0 

p3 pi sin (0 + $) cos (0 - <(>) 


where 


_ A 2 sin (t> cos 0 
s3 si sin (0 + <}>) 


For both the reflected and refracted waves, e is the angle of incidence and <$> is the angle 
of refraction obtained using Snell’s law. 

For normal incidence these equations can be written much more simply using the index 
of refraction, n. For the reflected wave: 

A _ » n - 1 
A 2 " A 1 n + 1 
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and for the refracted wave: 


A = A -2- 
A 3 ft l n + 1 

The attenuation (or absorption) of a transmitted wave can be found using the term, e 
where x is the distance traveled in the material with an attenuation factor of a. 

The above calculations for the electromagnetic response were coded into the CSTEM 
program. Another related feature that was developed is a method to determine what element 
faces of the model are exposed to an electromagnetic wave. This was done by determining how 
much of the edges of an element face can be seen from a viewpoint located on the -Z wave 
coordinate system axis looking in the direction of propagation. Once the exposed faces have 
been determined, an electromagnetic analysis can be performed for a wave impinging on each 
exposed face. 

The exposure analysis is done by considering the faces of an element where the tests are 
determined by the corner nodes only. The global coor dinates of the corner nodes are rotated 
to a local system so that the local Z axis lies parallel to the direction of propagation of the 
incoming electromagnetic waves. The face is split into two triangles and an area coordinate 
system is set up for each of the triangles. The exposed amount of each side of these triangles 
is determined by comparing with all other faces. The exposed amount of the face is assumed 
to be the same as the amount of the boundary of the face which is visible from the viewpoint 
of the incoming waves. 

The frontal area exposed to the waves can be determined by summing over the faces. The 
frontal area of each face is found as the area projection of the face on the local X-Y coordinate 
system plane (which is perpendicular to the wave propagation direction). This area is then 
modified according to the amount of the face visible to the incoming waves, and summed into 
the total exposed frontal area. 

The computer program WAVES has been successfully installed into the CSTEM analysis 
system. This program calculates the reflection and transmission of electromagnetic waves 
given a stackup sequence of materials and their electromagnetic properties. Similar to the 
other electromagnetic analyzers in CSTEM, WAVES will calculate a single cross section at a 
time. Multiple cross sections are handled as separate calculations as was done previously. The 
material properties needed are the real and complex parts of the relative permittivity and 
permeability. 

A description of the electromagnetic analysis portion of the CSTEM code as is currently 
implemented is included as Table 2. This description is a summary in words of the steps and 
the order in which they are followed in performing an electromagnetic analysis with the wave 
equations. The alternate table lookup procedure is not a part of this description. Included in 
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ROUTINE 


Table 2. Description of Absorption Code. 


ABMAN 

ABINP2 

ABMAN 

ABMAN 

ABMAN 

TRANIO 

ABMAN 

ABMAN 

(EXPOSE) 

ABMAN 

ELEMIO 

ABMAN 

ELGEQM 

ELDISP 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

S20IF 

JACOB 

ABSCAL 

ABSCAL 

LAraOT 

ABSCAL 

CROSS 

DOT 

ABSCAL 

CROSS 


CROSS 

DOT 

ABSCAL 

ABSTEM 

IMPEDE 

AB 

ABSCAL 

ABSCAL 

ABSCAL 

ABSCAL 

A10U8 

ABSCAL 

ABSCAL 

ABSCAL 

ABMAN 

1NEXT 

ELPATH 

ABMAN 

ABMAN 

ABMAN 

ABMAN 


Read waves and paths to be calculated (input sheets) 

Read # of absorption load cases to be done. 

LOOP over load cases. 

: Read wave -path combination for this load case 

; Get transformation defining wave coord system 
LOOP over paths 

: : Get starting element and face for this path 

(this is specified or from exposure analysis) 

: Get parameters describing incident wave 

(magnitude, freq, polar, angle, free space prop's) 

: Get element Information 

LOOP over elements in path 
; Get element geometry 

: Update nodal coords with displacements 

Calculate properties entering element 
: Check element face parallel to layers 

: : : Initialize thickness coordinate 

Set up unit vectors for propagation and E field 
directions, including polarization. 

LOOP over layers. 

Get layer mat'l properties, thickness 
: ; : Determine elem coord system axis thru thickness. 

: : (Integration loop here, but not really doing an 

integration. Calculations done at centroid.) 

: : : : Form shape functions (at centroid). 

: Compute Jacobian (relates elem system to global) 

Calculate global coordinates at layer upper and 
lower surface (at centroid) . 

: Calculate thickness of layer. 

: : Calculate global orientation of mat'l properties 

Check for same mat'l and orientation. If so, 
electric properties are same so attenuate only. 
Calculate an inward normal to impinged surface. 
Find angle between surface and impinging ray. 
Check for normal incidence (to save time). 

Form a normal to plane of Incidence. 

(plane of incidence Is defined by propagation 
direction and surface normal) 

: Set up an incidence coord system. 

Find angle between E and plane of incidence. 
Decompose E to in plane/out of plane components. 

: Interpolate electric properties for temp., freq. 

Calculate impedance (ETA). 

Calculate ph ase shift and attenuation constants 
(Integration loop ends here) 

Form index of refraction and angle of refraction 
Calculate amount of E reflected and refracted. 
Find distance travelled by refracted E (flat) 
Calculate amount refracted E is attenuated. 
Update thickness, exiting mat'l electric prop's. 
: : LOOP END (layers) 

Update outgoing E magnitude, orientation. 

Check for error (bad face or bad impinge angle) 

Find elements adjacent to present element. 

Find element connected to present element exit face 
LOOP END (elements) 

: Print out results for this path. 

LOOP END (paths) 

LOOP END (load cases) 
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the description are the subroutines where each step is performed. It is also easy to see the 
number of times a given step must be done in this type of analysis. 

2.3 CSTEM Acoustic Capability 

The capability was incorporated into CSTEM to determine acoustic characteristics due to 
structural vibration. The approach being to determine the radiation efficiencies of the 
structure for each vibration mode as a function of frequency. The structure geometry and mode 
shapes are obtained from an eignevalue run of a finite element model. Once the. radiation 
efficiencies have been found, the overall acoustic characteristics are determined by a 
summation of the contribution of the sound levels generated by each mode, given a prescribed 
force on the structure. Two methods for computing eigenvalues and eigenvectors have been 
implemented in CSTEM, the secant method determinant search and subspace iteration. These 
capabilities have been validated and verified against close form solutions. 

The CSTEM calculation of the sound power produced by a vibrating structure begins with 
the determination of the surface of the structure. The surface of the structure is found 
automatically from the element connectivities. In addition, a portion of the surface can be 
eliminated from a sound power calculations by masking it out. This masking is useful if only a 
portion of the surface is considered to be radiating sound. 

The eigenvalues and natural frequencies of the structure are first found. The radiation 
efficiencies for each mode are then found as a function of frequency based on the mode shapes 
and geometry of the structure. Once the radiation efficiencies are known, the sound power 
produced by a forced vibration of the structure can then be found by a modal summation of 
the contribution of each mode. 

The sound power calculated is based on the equation developed by Lord Rayleigh. This 
integral equation may be applied to extended plane surfaces. Other analysts have used this 
equation for nonplanar structures with success (Reference 1) for determining the radiation 
efficiencies and sound power. Radiation efficiency and sound power calculations used in 
CSTEM are based on References 2 and 3. 

The radiation efficiency calculations have been verified with a test case of rectangular 
plate. The results of running the model in CSTEM and that produced by Wallace are shown 

in Figure 22. 

2.4 CSTEM Tailoring 

The CSTEM tailoring capability is built on the STAEBL (Structural Tailoring Of Engine 
Blades) computer program from NASA Lewis. This was developed to perform engine fan and 
compressor blade numerical optimizations. 

The program consists of two major modules; CONMIN, which performs the optimization 
and ANALIZ, which supplies values for the parameter involved in the optimization. 
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Sound Power Calculation Input Description. 


UNE 

NUMOMG 

NDRMSK NODDIR NODFRC 

NUMOMG 


Number of frequencies that radiation efficiencies and sound power are 
calculated 

NDRMSK 


Number of mask points 

NODDIR 


Direction of nodal forcing function 
X direction = 1 
Y direction = 2 
Z direction = 3 

NODFRC 


Node number at which force is applied 

UNE 

FROMIN 

FROMAX 

FROMIN 


Minimum frequency for sound power calculation (Hz) 

FROMAX 


Maximum frequency 'or sound power calculation (Hz) 

The NUMONG frequencies will be logarithmically distributed between 
FROMIN and FROMAX. 

UNE 

FORCE 


FORCE 


The magnitude of the applied force (Ibf) 

UNE 

ETA(1) 

ETA(2) ETA(NEIG) 

ETA(i) 


The modal loss factor for each mode. One value required for each 
eigenvalue calculated. 

UNE 

DIRMSK(I) 

DIRMSK(2) DIRMSK(3) 

DIRMSK(i) 


The X, Y, and Z coordinates of the masking point. 


A masking point is used to eliminate some of the surfaces of elements 
from being included in the sound power calculations. Faces of 
elements interior to a structure are automatically eliminated by the 
program and need not be considered for masking. 

It Is often desirable to not include all surfaces of a structure for sound 
calculations. For example the exterior surface of a pipe might be the 
only surface desired for sound power calculations. The interior of the 
pipe can be masked out by specifying a masking point anywhere in 
the interior of the pipe. 

Masked surfaces are defined to be any element surface which has a 
positive dot product between the surface of an element and a vector 
from the element face to the masked point. Thus only one point is 
required to mask all interior surface faces of a straight pipe. 

Output of radiation efficiencies are dimensionless and the sound power is output In WATTS. 
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EFFICIENCIES 
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Figure 22. Correlation Between CSTEM Predictions and Closed Form Solution for the 
Radiation Efficiencies of a Rectangular Plate. 





CONMIN works on a set of up to 125 global variables to perform the optimization. A 
subset of these are selected by input and designated as design parameters, constraints, or object 
functions. These variables are quantities such as blade thickness, layer orientation, root chord 
length, and blade weight. 

A main driver program, COPES, calls CONMIN initially to set up the global variables. It 
then calls ANALIZ which reads input and computes a sufficient set of parameters for 
optimization to precede. CONMIN supplies ANALIZ with a modified set of input and the 
procedure continues until some convergence criteria is met. 

The CONMIN module was abstracted from STAEBL and coupled with CSTEM 
structural, thermal, electromagnetic, and acoustic application modules. This was successfully 
debugged and the following verification/validation cases have been run. 

The first case involved the tailoring of the natural frequencies of a plate. In this, seven 
global variables were used. Three of these were the first three natural frequencies and the 
other four were ply angles of the material that made up the model. The model used was a 
cantilevered plate consisting of nine 20-noded elements with a total of 96 nodes. The material 
was a 4 ply composite material. 

In the optimization input, the lowest frequency was taken as the objective function to be 
minimized by varying the 4 ply angles and constrained by the other 2 frequencies. After 10 
iterations, the object function was reduced from 565 CPS to461 CPS. Each iteration took about 
10 seconds on a Cray XMP. 

1 

The CSTEM tailoring program is capable of tailoring a finite element model to control 
the noise produced by the structure. 

A 96-noded 9-element square plate model was used as a test case. The noise was in 
response to 4000 DB 565 cps source. The plate was constructed of 4 plies of composite material 
at different orientation angles. These angles were taken as design variables. Their variation 
leads to the change to natural frequencies and associated mode shapes of the plate. 

As shown in Table 3 the noise level (the objective function) was reduced by a factor of 33 
after 24 iterations (200 seconds on the CRAY XMP). 


Table 3. Sound Power Tailoring. 


Iteration 

Number 

Ply Angles 

Frequencies, cps 

Noise, watts 

1 

45.0 

90 

90 

45.0 

5650 

1447 

3074 

7.1 x 10j 


67/5 

90 

90 

45 

674 

1403 

2892 

1.8 x lO" 6 

1 

45 

90 

90 

67.5 

674 

1403 

2892 

1.8 x 10"! 


67.5 

45 

90 

67.5 

734 

1435 

2782 

4.8 x 10' 7 

15 

39 

90 

90 

79 

734 

1256 

2721 

2.2 x 10"® 

20 

79 

90 

90 

39.6 

734 

1256 

2721 

2.2 x 10" 6 

23 

90 

90 

90 

60 

1007 

1377 

2437 

2.16 x 10- 7 

24 

79 

90 

90 

79 

981 

1393 

2473 

2.13 x 10‘ 7 
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An example tailoring problem was run with the CSTEM electromagnetic analyzer using 
commercially available electrical properties for composite materials (see Table 4). These 
properties are for a single frequency and temperature only and are given in a different form 
than needed by the electromagnetic analyzer. It is assumed that these properties are typical of 
those that the CSTEM analyzer may be used for, so the program was modified to accept 
electrical material data given as dielectric constant and loss tangent. The permittivity, 
permeability, and conductivity needed by the electromagnetic analyzer can be obtained from 
these properties. 

Tabic 4. Electrical Properties Comparison. 


Material 

Fibar/Raain 

Lota 

Dielectric 

Constant 

Tangent, 

Spectra/Polyester 

2.20 

7 

Spectra/Epoxy 

2.12 

20 

Spectra/Vinylester 

2.16 

21 

Spectra/Polyethylene 

1.99 

17 

Glass/Polyaster 

4.30 

150 

AramkJ/Polyester 

3.50 

500 


The permeability when not mentioned is assumed to be that of free space (1.257E-6 
henry/meter or 3.192E-8 henry/inch). The permittivity can be found using 


e = Er Co 

where t r is the given relative dielectric constant and eo is the permittivity of free space 
(8.854E-12 farad/meter or 2.249E-13 farad/inch). The conductivity is found from the equation 

for the loss tangent as 



where <o is the frequency and a is the conductivity. 

The amount of electromagnetic energy transmitted through a 1/2 inch thickness of material 
was calculated for a series of impingement angles. The results arc shown in Figure 23. It can 
be seen that the greatest amount of transmitted energy is at normal incidence. The difference 
in material types is also evident. 

Table 5 gives the derived electromagnetic properties e and o(p « 1.0) for these materials. 
The range in properties, particularly since a is a function of frequency, will make for significant 
optimization problems. 
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5 INCH THICKNESS 



Angle From Normal Incidence, Degree 
Q SPPUtTH + SPEPOXY O OLPOLY 

Figure 23. Amount of Electromagnetic Energy Transmitt* J. 


Table 5. Derived Electromagnetic Properties. 


Material 

Fiber/Retin 

Permittivity 
e(1 0 1 *) 
farad/inch 

Conductivity 
at 10 GHz 

o(10' 7 ) Siemens/inch 

Spectra/Polyester 

4.95 

34.36 

Spectra/Epoxy 

4.77 

95.36 

Spectra/Vinyiester 

4.86 

102.01 

Spectra/Polyethylene 

4.48 

76.08 

Glass/Polyester 

9.67 

1450.61 

Aramid/Poiyester 

7.87 

3935.75 

- 


As a test case for electromagnetic absorption, a cantilever beam of eight 8-noded bricks 
was used. The beam was constructed of two layers of different materials (orthogonal to the 
X-Y plane). These materials were a SPECTRA fiber in an epoxy matrix (SPEPOXY) and 
AJRAMID fibers in a polyester matrix (ARAPOLY). Two design variables (the percentage of 
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thickness of each material) were used. The absorption in the Z direction was the objective 
function to be maximized. The results are as shown in Table 6. 


Table 6. Electromagnetic Absorption Tailoring. 


SPEPOXY, % 

50 

54 

45 

25 

33 

ARAPOLY. % 

50 

46 

55 

75 

77 

Absorption, % 

9.3 

8.63 

10.14 

13.48 

13.82 


In an effort to determine how to improve the efficiency of the CSTEM optimizer, 
discussions were held with Dr. J.K. Casey, an applied mathematician at GE who has experience 
working in optimization problems. One concern was that a great deal of comp utational effort 
would be required in order to optimize all of the different aspects of a CSTEM problem. 

Dr. Casey pointed out that much of the efficiency in an optimization problem in the hands 
of the user of the program. To efficiently optimize a problem the user should decompose it 
and optimize a single variable separately in several one dimensional problems rather than try 
to optimize several variables in one multiple dimension problem. The order in which the 
optimization is carried out can be chosen to maximize the performance gained when 
optimizing a variable at a time. The user should choose to first optimize those variables that 
experience shows to be most critical. 

Also, a simple, coarse mesh should be used when optimizing. If a problem is sensitive to 
mesh density it is possible to use a coarse model to get derivative values that the optimizer 
needs and a finer model to compute the results. 

Presently, the optimizer in CSTEM can only work on a single objec t function. This is the 
same as the STAEBL program that the optimization package in CSTEM comes from. A way 
to optimize multiple object functions is being examined. This would require an additional 
executive program around the optimization loop. 

One possibility would be to direct the optimizer to perform an iteration on each of the 
object functions and store results from the separate optimizing iterations. An additional 
decision module would also be required to determine the averaged direction to take which 
would make some improvement in all of the object functions. This would continue until 
arriving at a point which requires opposite changes in order to improve the different object 
functions. 

Another way to look at multiple object functions might be to optimize each separately, 
using previously optimized object functions as constraints. This would require prioritizing the 
object functions, iterating on the most important one first then using it as a constraining value 
while iterating on a second object function. It should be pointed out that our meeting with Dr. 
Casey indicated that it may be better to decompose the problem and optimize each objective 
as a separate problem. 
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To prepare for more complicated optimization work a large 8-noded brick model of a J-85 
compressor blade was recactivated (Figures 24 and 25). This model consists of 100 8-noded 
isoparametric elements and 242 nodes. The bill-of-material for this blade is titanium. To 
benchmark the model and further verify the large deformation logic in CSTEM, comparison 
runs were made among the centrifugal stiffening predictions of NASTRAN (MSC), MASS 
(GEAE), and the updated Lagrangian approach in CSTEM. A rotor speed of 16,000 rpm was 
assumed and titanium material properties. The lowest three natural frequencies were 
computed by each program as shown in Table 7. 

Table 7. Correlation of CSTEM Nonlinear Eigenanalysis. 


Mode 

CSTEM 

Frequency, cps 

MASS 

Frequency 

NASTRAN 
Frequency, cps 

1 

936 

935 

981 

2 

2763 

2763 

2092 

3 

3265 

3278 

3302 


CSTEM correlates excellently with MASS which has been related with test data 
NASTRAN gives consistently stiffer predictions. 

The use of CSTEM in the analysis of these fairly large blade mode’s i.as pointed out the 
need for an alternative method to generate element layers in a model. The present manual 
layer input is an element by element specification which requires the user to manually 
determine where element boundaries lie within a layup. This can be difficult. 

Another method has now been included in CSTEM which allows the usci to specify a layup 
and then apply it to a cross section of elements. To specify the cross section the user need only 
input the surface element and the element axis that points through the thickness (along the 
cross section). The program will determine what elements are in the cross section and the 
thickness of the elements at their centroids. Using the element thicknesses and the thicknesses 
of the layers comprising the Lyup the program determines how many layers lie within each 
element. These element layers may then be copied to other similar cross sections by specifying 
the beginning element of each of these cross sections. 

A UDF blade model has been acquired for CSTEM development. This propulsor blade 
model consisted of 756 nodes and 480 eight-noded isoparametric elements. The model was 
run previously on our in-house program MASS using orthotropic material properties. This 
model was run on the CRAY version of CSTEM. The composite properties were simulated 
by 23 different materials using 481 reference nodes to define the composite directions. The 
large displacement analysis capability of CSTEM was used to account for centrifugal stiffening. 
Six load increments were used to spool up from flight idle to cruise. At the cruise condition, a 
deformed position eigenanalysis was performed using the subspace iteration technique. The 
lowest three natural frequencies were obtained and are compared in Table 8 against MASS. 
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Table 8. De formed Position Eigenanalysis 
With Subspace Iteration. 



CSTEM 

Mass 

1 . 

281 

278 

2. 

593 

598 

3. 

858 

853 


This required approximate 370 CPU seconds on the CRAY. 

2.5 CSTEM Composite Micromechanics Using ICAN 

Composite micromechanics analyses in the CSTEM program are performed using a 
modified version of the ICAN computer program obtained from NASA-Lewis. This program 
and the inputs are described in the ICAN User’s Manual. In the CSTEM application these 
inputs are obtained from the CSTEM input deck, are calculated from the results of the CSTEM 
finite element program, or are defaulted to a value such that a certain type of analysis is always 
performed. The following is an attempt to describe the ICAN installation in the CSTEM 
program and how the specific ICAN inputs are obtained or calculated. 

The ICAN routines can be accessed in two separate ways. This results from the fact that 
the ICAN program itself incorporates routines from a program called INHYD which calculates 
layer properties from the constituent fiber and matrix properties. The" * 2AN/INHYD type 
routines can be accessed to generate elastic material properties for composite systems whose 
constituent properties are contained in a data bank. Section X.4 in the CSTEM input sheets 
describe the input variables for this type of use. The actual micromechanics analysis can be 
performed as well using the ICAN routines. In this usage, the elastic material properties must 
be obtained through the use of the ICAN/INHYD routines contained in I--. >N as mentioned 
above. Since ICAN performs its analysis for conditions at a cross section of the laminate, 
additional input specifying this location is needed. Section XVIII. in the CSTEM input sheets 
describe the inputs specifying these cross sections. 

The inputs needed to use .CAN to generate elastic properties are those contained in the 
Ply Details cards and Material System Details cards described in the ICAN User’s Manual. 
The material identification number is assigned consecutively in order of input. All properties 
are returned in the local material coordinate system and are incorporated into the global 
system by the CSTEM program so that the orientation angle and thickness in the Ply Details 
card group are not necessary. The use temperature is passed from the CSTEM program to the 
ICAN/INHYD routines and is interpolated from nodal temperature values. There is a check 
for a trash temperature value of 1234567. in which case a reference temperature of 70°F is 
used. At present this would occur only if using the linear heat transfer analysis option without 
a corresponding structural analysis. 

The ICAN materials data bank is assumed to exist as a file called ’ICANBNK’ on the 
machine where the CSTEM program is being run. This file is opened by the CSTEM program 
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(subroutine RDCAN) when the option for ICAN generation of material properties is 
activated. 

The use of ICAN as part of the CSTEM program for micromechanics analysis requires 
only that the material properties be generated by ICAN as previously described, that the 
location(s) where the analysis applies be specified, and whether strains and curvatures or loads 
are to be calculated from the finite element results as input for the ICAN analyzer. The ICAN 
analysis is done immediately after the structural results have been computed and printed out 
so that an ICAN load case corresponds to a structural load case. However, an ICAN analysis 
can be done for up to 10 different cross sections in a load case. 

The cross sections within the structure where the ICAN analysis is to be performed are 
specified by listing the element numbers through the cross section starting from the bottom. 
The calculations to obtain the necessary information for ICAN are done at the centroids of 
the layer surfaces within the listed elements. Therefore, the cross section element centroids 
should be lined up normal to the layers. A maximum of 10 elements are allowed in defining a 

cross section. 

The n um ber of layers in the cross section are determined by counting the number of layers 
in each consecutive element of the cross section starting from the bottom of the laminate (the 
first element listed). There is no check for layers extending between elements, so if a layer lies 
in two adjacent elements it would be treated as two layers having the same orientation. The 
present limit on number of layers through the cross section is 100. 

The limits on number of cross sections, elements through the cross section, and layers 
through the cross section are all contained in PARAMETER statements. Therefore, they can 
be fairly easily changed by modifying this PARAMETER statement and recompiling the 
CSTEM program. 

The analysis options set by the booleans described in the ICAN User s Manual are 
hardwired in the program (Table 9). 


Table 9. ICAN Boolean Settings. 


COMSAT 

True 

BIDE 

True 

CSANB 

False 

NONUDF 

True 


The boolean RINDV is determined by the CSTEM input variable ICAN. RINDV is True 
if ICAN = 1 and False if ICAN = 2. 

Much of the data required by ICAN as described in the Ply Details cards and Material 
System Details cards of the ICAN User’s Manual is already available from the generation of 
the material properties using the ICAN/INHYD routines. The percent of moisture of each 
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layer, which was input for the generation of material properties, is transferred to storage in 
array PL, row 72 for use by ICAN. 

The use temperature for each layer is interpolated from the temperature of the element 
nodes to the midsurface of the layer using the element shape functions. ICAN uses the 
difference in this use temperature and the cure temperature of the material. This value is 
stored for each ply in array PL, row 50. 

The thickness of each layer is calculated by first locating the centroid of the upper and 
lower surfaces of the layer using the element shape functions and nodal coordinates then 
calculating the distance between these two points. This value is stored for each layer in array 
PL, row 7. 3 

The orientation angle of the layer is obtained from the element layer information stored 
using the LAYIO routine of the CSTEM program. This angle is placed in the THLC vector 
used by the ICAN routines. 

The bulk of the calculations done using the finite element results are to obtain the loadings 
for the ICAN analyzer. These calculations assume that the strain varies linearly through the 
cross section, while the stress may be discontinuous. The loadings can be either strains and 
curvatures or loads in the form of stress resultants and couples. Strains and stresses at the cross 
section element integration points, along with the element geometry information are the 
quantities used from the finite element analysis. The integration points are generally in a 
Gaussian distribution on each layer midsurface. In the case of an element with only one layer 
the integration points will be in a three dimensional Gaussian distribution. 

For the strain and curvature loading option of ICAN, the strains and curvatures at the cross 
section midsurface (reference plane) are needed. The strains are calculated by first calculating 
the thickness through the entire cross section as the distance between the centroid on the upper 
surface of the cross section and the centroid on the lower surface. These centroids are 
calculated using shape functions and nodal coordinates. Once the cross section thickness has 
been determined the layer thicknesses are summed beginning from the bottom surface until 
this sum exceeds half of the cross section thickness. The strains at the integration points of the 
element in which the reference plane lies are then interpolated to the reference plane centroid 
using Lagrange interpolation functions. 

The reference plane curvature is calculated using the equation for strain from laminated 
plate theory which can be written as: 

{e} = {e°}-Z{K°} 

where {e°} is the reference plane strain, {K°} is the reference plane curvature, and Z is the 
distance from the reference plane to the location where the strain, {e}, is desired. This equation 
can be solved for the curvatures. 

{K°} = l/Z({e°}-{e}) 
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The curvatures can be calculated using this equation, the strain at the centroid of the fop 
layer of the cross section, aid the distance from the reference plane to this centroid. 

The in-plane stress resultants or membrane loads used by ICAN can be calculated using 
the equations 


M2 

N x = o x dz 

J-h/2 

Ml 

N v = Gy dz 

' J-h/2 

j.h/2 

dz 


while the equations for the stress couples or bending resultants are 



j.h/2 

- 

M x = C x z dz 

x J-h/2 

- 

j.h/2 


M v = Gy z dz 

J J-h/2 - 

: 

j.h/2 • 


M *> - J _ W2 a *y z dz 


and the equations for the transverse shear resultants are 
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r h/2 

Ox — I °xz ^ 2 . 

J- h/2 

j.h/2 

Qy I ^yz ^ 

Jh/2 7 

In these equations, h is the cross section thickness and Z is the distance from the reference 
plane. 

Since the stress through the cross section is assumed to be piecewise linear these equations 
must be integrated piecewise, integrating over each individual layer and summing the layer 
results over the cross section. For the membrane loads this can be written in general as 


layers 


f ' + 7 

{o} dz 


• 4 - 


2 


where Zr is the distance from the reference plane to the layer midpiane and h is the layer 
thickness. Assuming the stress to be linear over the layer it pan be written 

= ^mid) + |f ((0 t0p } - {O^}) 


where {amid} is the layer midsurface stress, {mop} is the stress at the layer top surface, {obot} 
is the stress at the layer bottom surface, and Z is the distance from the layer midsurface to a 
point in the layer and can be written 

Z1 = Z - Z, 


Substituting this stress function into the integral gives 
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layersl 



he 

rZr+ 2 

, t a iop ~ a boJ 

dz + t-r zdz - 

h e 

^^top ^bot^ 

he ^ 

hi 1 
r z r+ y 

dz 


he 

z - — 


K 


r 2 

” J 


39 



wiich when integrated results in the following general equation for the membrane loads. 

{N} = X h > 

layers 


The transverse shear resultant integrals are of the same form as those for the membrane 
load^e ^ndmg resukits are somewhat different and use of the same procedure results 

in a general equation for the bending resultants. 


(®top ^botl 


X , 1. '-'top 

{OmidlZph^- ^ 




layers 


th, n., entities needed in the equations for the loads are obtained by looping through all 

dielayetsbegi^ngfrom^^bo^^o^t^ CTOsssectwi^calCTilmingth^®m^MfiQn^oftha^ 

S5Ha£^i#fflSM=s 

so that the stress at These points can then be obtained from the strains. 

2.6 Multiple Layer Elements In CSTEM 

The information necessary to define the layering of an element is the thiclmess axis, the 
layer material, the layer thickness, the number of layers, and the orientation oft y 

mmmmm 

toes he at *e“ 1 Sd 1 locations of its perpendicular rst axis. The rst system defiruuon for 
CSTEM is shown in Table 10. 
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Axis 

Comer Nodes 

+t 

1 

2 

3 

4 

-t 

5 

6 

7 

8 

+ s 

1 

2 

6 

5 

-s 

3 

4 

8 

7 

f T 

1 

4 

8 

5 

-r 

2 

3 

7 

6 


* 

so that going through the element thickness means we are going perpendicular to the layers. 
This way of defining element layers means that when modeling a layered structure the mesh 
must follow the layer surfaces so that layers run between opposite element faces. Layers can 
not cut diagonally across an element. 

The layer material is defined by a material number and the properties associated with this 
material are stored by reference to this material number. Layer materials are treated as 
orthotropic materials, so isotropic materials can be used since they are a subset of orthotropics. 
Orthotropic materials have directionality so a material coordinate system is associated with 
each material. The material properties are specified with respect to this material coordinate 
system with the material X or 1 axis usually being the principle or stronger direction, although 
this is not a requirement. 

Since the layers are of constant thickness in isoparametric space, the layer thicknesses are 
stored as fractions of the element thickness. Therefore, if the element varies in thickness then 
so do the layers. CSTEM is set up so that layer number 1 is that the negative end of the thickness 
axis. The layers are then numbered consecutively progressing along the tnickness axis in the 
positive direction. 

There are two ways to define the orientation of the layer with respect to the global 
coordinate system. One is to define a transformation directly from the material coordinate 
system to the global coordina '.e system. In CSTEM a particular transformation has a number 
associated with it much like the materials. For global orientation of element layers the 
transformation number is all that is required and the transformation associated with that 
number is used to rotate the material properties directly to global. 

Many structures are of complex shapes and curvatures which makes it difficult to specify 
the material orientation with respect to a global coordinate system. Since the elements in a 
mesh generally follow the shape of a structure (as do the layers) it may be much easier to specify 
the material orientation with respect to the element coordinate system. The material 
orientation with respect to the global system can then be calculated since the information 
relating the rst system to the global system is available from the element geometry definition. 
Using this method in CSTEM assumes that the material coordinate system is aligned so that 
the material Z or 3 axis is coincident with the element thickness axis. For local orientation of 
element layers the angle of rotation about the material 3 axis is all that is required. The 
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reference axis for this rotation angle depends on which of the rst axes is defined as the thickness 
axis. Table 11 shows the reference (0 degree) orientation of the material coordinate system to 
the element rst system for each thickness axis case. 


Table 11. Unrotated Material Orientation with Element 



Material 

Element 


Axle 

Axes 


1 

rst 


2 

s t r 

(thickness axis) 

3 

t r s 


The orientation of the element coordinate system with respect to the global coordinate 
system is contained in the Jacobian matrix. The Jacobian contains both the angular relation of 
the element coordinate system to the global coordinate system as well as the stretching or 
shrinking of the element axes in relation to the global axes. For the material orientation only 
the angular relationship is required. Therefore, the Jacobian must be normalized row by row 
to get the transformation matrix from the element coordinate system to the global coordinate 
system. (If the Jacobian were to be normalized column by column the transformation matrix 
from the global coordinate system to the element coordinate system would be obtained). A 
possible drawback in using the Jacobian in this way is that the rst coordinate system may not 
be an orthogonal system when viewed in the global coordinate system. This would occur when 
element comers are not right angles. 

The rotation of the material matrix from the material coordinate system to the global 
coordinate system is calculated as: 

id ] 8 = m T pi 1 m 

Therefore, to get the complete transformation matrix from the material coordinate system 
to the global coordinate system the normalized Jacobian (which rotates from element to 
global) is premultiplied by the planar trans forma tion matrix which rotates the material 
coordinate system to the element coordinate system. 

Further example problems using composite materials are being run with the structural 
portion of the CSTEM code to gain experience with the layered element. In many composite 
structures the number of plies through the thickness varies along the structure/Test cases are 
being prepared to determine the best way to model this ply dropoff using the CSTEM layered 
element, which must contain the same number of layers throughout the element. 

Figure 26 shows an example of a cross section geometry which changes from four layers 
through the thickness to two layers through the thickness. Two m ethods can be used to mod el 
this with multiple layer elements. In the first method, the elements are tapered to follow the 
geometry of the structure. This results in a mismatch of layers between one element to another 
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(Figure 27). This poses no computational difficulties since the stiffness for each element is 
generated separately. In the second method the element geometry would not follow the actual 
structure geometry, but would be modeled as if there was no drop off of plies. The ply drop 
off is modeled by zeroing out the material properties of the absent layers (Figure 28). Again 
this does not cause a computational problem since the layers that are present will add their 
stiffness to the element. A combination of these two methods might be the most effective in 
a general problem. 

A tapered composite cantilever beam example was used initially to examine this problem. 
A model with only one layer per element was used to validate the results of modeling this beam 
with multiple layer elements (see Figures 29 and 30). 
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Figure 26. Actual Cross Section Geometry. 
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Figure 23. Zero Material. 
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3.0 Conclusions 

The first five tasks of the NASA Statement of Work has been successfully completed. 
Unique graded material finite elements have been developed with structural, thermal, 
electromagnetic, and acoustic capabilities. A CSTEM analyzer has been developed which 
performs coupled structural/thermal/electromagnetic/acoustic analysis. Other advanced 
capabilities are the large deformation analysis of graded composite structures, deformed 
position eigenanalysis, transient nonlinear heat transfer with right- hand-side pseudo-fluxes, 
electromagnetic absorption, and sound power predictions. A CSTEM tailoring capability has 
been developed by marrying the optimization routines from the NASA supplied STAEBL 
program with the CSTEM application routines, Verification and validation cases have been 
run. One final task remains, the CSTEM tailoring of simulated components. 
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